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LEAST SQUARES SHADOWING METHOD FOR SENSITIVITY 
ANALYSIS OF DIFFERENTIAL EQUATIONS 

MARIO CHATER*, ANGXIU NI*, PATRICK J. BLONIGAN*, AND QIQI WANG* 

Abstract. For a parameterized hyperbolic system ^ = f{u,s) the derivative of the ergodic 

average (J) = lim 7 ’_j.oo ^ s) to the parameter s can be computed via the Least Squares 

Shadowing algorithm (LSS). We assume that the sytem is ergodic which means that (J) depends 
only on s (not on the initial condition of the hyperbolic system). After discretizing this continuous 
system using a fixed timestep, the algorithm solves a constrained least squares problem and, from the 
solution to this problem, computes the desired derivative . The purpose of this paper is to prove 
that the value given by the LSS algorithm approaches the exact derivative when the discretization 
timestep goes to 0 and the timespan used to formulate the least squares problem grows to infinity. 
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1. Introduction. Consider the differential equation parameterized by 5 C M 
and governing u{t) G : 


/ f = /K^) 

\ u(0) = Uq Uq € 


( 1 . 1 ) 


The differential equation is assumed to be uniformly hyperbolic (details in section 3). 
We are also given a cost function J(u,s) : x R —>■ R and assume that the 

system is ergodic^ i.e., the infinite time average : 

(-^)('S) = 4/ J{u{t),s)dt (1.2) 

T->+oo T Jq 

depends on s but does not depend on the initial condition m( 0). The differentiability 
of (J) with respect to s has been proven by Ruelle [1]. Obtaining an estimation of 
is crucial in many computational and engineering problems. Indeed, many ap¬ 
plications involve simulations of nonlinear dynamical systems that exhibit a chaotic 
behavior. For instance, chaos can be encountered in the following fields : climate and 
weather prediction [2], turbulent combustion simulation [3], nuclear reactor physics 
[4], plasma dynamics in fusion [5] and multi-body problems [6]. The quantities of in¬ 
terest are often time averages or expected values of some cost function J. Estimating 
the derivative of (J) is particularly valuable in: 


• Numerical optimization. The derivative of (J) with respect to a design pa¬ 
rameter s is used by gradient-based algorithms in order to efficiently optimize 
the design parameters in high dimensional design spaces (see [7]). 

• Uncertainty quantification. The derivative of (J) with respect to a pa¬ 
rameter s gives a useful information for assessing the error and uncertainty 
in the computed {J) (see [8]). 
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For example, we could obtain a useful information about the impact of mankind on 
the climate by computing the derivative of the long time averaged global mean tem¬ 
perature to the amount of anthropogenic emissions ([9] shows how sensitivity analysis 
is used in climate studies). In the simulation of a turbulent airflow over an aircraft, 
estimating the derivative of the long time averaged drag to a shape design parameter 
is of extreme importance for engineers allowing them to improve their design [10]. It 
has been shown that in many of these practical examples, the quantities of interest 
exhibit ergodic properties, popularly known as chaotic hypothesis [11], [12]. 

When it comes to computing , conventional methods based on linearizing the 
initial value problem (1.1) become ill-conditioned when the system is chaotic. They 
compute derivatives that are orders of magnitude too large and the error grows ex¬ 
ponentially larger as the simulation runs longer [13],[14]. This failure is due to the 
so-called butterfly effect and the explanation has been published by Lea et al. [13]. 
Some algorithms have been developed to overcome this failure. Lea et al. proposed 
the ensemble adjoint method which applies the adjoint method to many random 
trajectories, then averages the computed derivatives [13], [15]. However, the algo¬ 
rithm is computationally expensive even for small dynamical system such as Lorenz’s 
one. Based on the fluctuation dissipation theorem, Abramov and Majda provided an 
algorithm that successfully computes the desired derivative [16]. Nonetheless, this al¬ 
gorithm assumes the dynamical system to have an equilibrium distribution similar to 
the Gaussian distribution, an assumption often violated in very dissipative systems. 
Recent work by Cooper and Haynes has alliviated this limitation by introducing a 
nonparametric method to estimate the equilibrium distribution [17]. More methods 
have been developed to compute in particular the Least Squares Shadowing 

(LSS) algorithm which computes it by solving a constrained least squares problem 
[14]. The big advantage of this method is its simplicity since the least squares problem 
can easily be formulated and efficiently solved as a linear system. Compared to the 
previously presented methods, LSS is less sensitive to the dimension of the dynamical 
system and doesn’t require any explicit knowledge about its steady-state distribution 
in phase space. 

This paper provides a theoretical foundation for LSS by proving that it gives a 
useful estimation of when the dynamical system is a uniformly hyperbolic flow. 
Compared to the discrete case (uniformly hyperbolic map) for which we already have 
a proof of convergence [18], the continuous case is more difficult to deal with due 
to the apparition of the neutral subspace (details in section 3). However, it is very 
important to treat the continuous case since most applications and real-life problems 
require a continuous description of the physics and involve differential equations. 

In the next section, the mathematical formulation of convergence is introduced 
as well as theorem LSS which will be proved in the following sections. Section 3 
presents the concept of uniform hyperbolicity for the readers who are not familiar 
with the subject. Section 4 points out the new behaviour and properties that come 
with continuous dynamical systems (in opposition to discrete maps). Section 5 defines 
the shadowing direction and proves its existence as well as uniqueness. Section 6 shows 
that the derivative of (J) can be computed using the shadowing direction and bounds 
the upper error. Section 7 then demonstrates that the least squares problem gives a 
good approximation of the shadowing direction. Finally, section 8 uses all the previous 
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results and concludes the proof of theorem LSS by showing that the estimation error 
vanishes as the least squares problem increases in size. 

2. Discretizating the problem and LSS theorem. In order to obtain an 
algorithm of practical relevance, a discrete version of the above problem should be 
formulated. First, we replace the differential equation (1.1) parametrized by s by a 
family of operators : let ips {u,h) : R™ x R+* —>■ R™ be the family of maps parametrized 
by s G R such that it is a "discretization" of the differential equation using a uniform 
time stepsize of h. In other words, if {u(t),t G {—oo, +oo)} is a trajectory satisfying 
the initial differential equation (1.1) for a particular s G R, we have : 

ips{u{t),h) = u{t + h) Vt G M 

In this case, (fs{',h) corresponds to a perfect numerical integration scheme with a 
stepsize of h. We ask for the differential equation to be "smooth" enough so that the 
maps ■) are C^. 

Then, assuming that all trajectories {u(t),t G R} belong to a compact set A and 
that s also lies in a compact set ^ C R, we can approximate ^ J(u(t), s)dt with a 

Riemann sum and have the following bound on the integration error: 


1 

T 


: h I 


J(u(t), s)dt— > 


i=0 


ih),s)) < /isnp(||HJ(-,s)||^)snp(||/(-,s)||^) 
seS seS 


( 2 . 1 ) 

where (ps{u{0), 0) = u(0) and ||/(•, s)||^, \\DJ{-, s)||“ are the infinite norms of /(•,«) 
and the derivative of T(-,s) with respect to the first variable on the compact set A 
respectively. Since the bonnd doesn’t depend on T, we finally have : 


-'^hJ{(ps{u{0),ih),s) = 0{h) 

T^oo 1 


2 = 0 


( 2 . 2 ) 


To simplify the expressions, we introdnce the new notation (ps(u(0),ih) = or 

even (ps{u{0),ih) = = uf = Ui depending on which parameter is fixed and 

when there is no ambignity. 


For a sequence {ui, i = 1, - ■ ■ , [^]} satisfying Ui+i = <Ps{ui, h), the Least Squares 
Shadowing method attempts to compute the derivative via 

Theorem 2.1 (THEOREM LSS). Under ergodicity and hyperholicity assump¬ 
tions, 


d{J) 

ds 


(s) = 1™ 1™ \iDJiui,s))vl^’'^^ + dsJiui,s) + {r]j’"'^\j{u^,s) - (J)(s))) 

Ai—>-0 1 —¥oo ± \ ^ L 


2=1 


^ [{DJ{ui,s))vl^’'^^ + dsJiui,s) + {r]j’"'^\j{ui,s) - (J)(s))) 


2=1 
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where £ R'” x R, i = 1,..., [^] is the solution to the constrained least 

squares problem : 

i=i 

s.t. = {Dips{ui,h))vj'^’'^^ +ds(ps{ui,h) + dhips{ui,h), 

where a is any positive constant and 11.11 is the Euclidean norm in R™. 

Here the linearized operators are defined as : 


{DJ{u, s))v 
{Dtfsiu, h))v 


(DyJ)(u, s) := lim 
£->■0 


J{u + ev, s) — J(u, s) 


{Dv(ps){u, h) ■= Im 


(ps{u + ev, h) - ipsiu, h) 


, , J{u,s + e)- J{u,s) 

OsJiu, s) := lim- 

e^O e 

a f r V’s+e{u,h) - ips{u,h) 

ds^Psiu, h) := hm- 

e ->0 e 

a f h\ r f8{u,h + e)-ips{u,h) 

dhips{u, h) := hm- = f{ips{u, h)) 

£->■0 e 


(2.4) 


{D J),{dsJ),{Dips),(dsips) and {dh^Ps) ^^re a 1 x m vector, a scalar, an m x m matrix, 
an m X 1 vector and an m x 1 vector, respectively, representing the partial derivatives. 


3. Uniform hyperbolicity. Let us now proceed to the presentation of the uni¬ 
form hyperbolicity property : we say that the dynamical system (1.1) has a compact, 
global, uniformly hyperbolic attractor A C R™ at s if for all t the map Psi', t) satisfies: 

1. For all uq £ R"*, dist(A, zi(t)) 0 where dist is the Euclidean distance in 

R™. 

2. There is a C £ (0,-l-oo) and A £ (0,1), such that for all u £ A, there is a 
splitting of R™ representing the space of perturbations around u : 

R"* =U+(u)0U-(u)0U°(u) (3.1) 

where the subspaces are : 

• V~^{u) := {u £ R"*/ \\{D(ps{u,t))v\\ < CA“*||u||,Vt < 0} is the unsta¬ 

ble subspace at u, 

• V~{u) := {u £ R™/ \\{D(ps{u,t))v\\ < UA^IIuHjVt > 0} is the stable 

subspace at u. 

• V°{u) := {af{u,s),ya £ R} = {adh<Ps{ui-i, h),ya £ R} is the neutral 
subspace at u. 

V~{u),V'^{u) and U°(u) are all continuous with respect to u. 

li r = r'^ -\- r~ -\- r^ with r’*' £ V~^{u), r~ £ V~(u), £ V^{u) and u £ A, the 

continuity of the three subspaces and the compacity of A implies that: 

. ||r+ 0 r“ 0 r°|| 

u,r+“-,r-o max(||r+||, ||r-||, ||r0||) 


/3 > 0 


(3.2) 
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This is because if /3 = 0, then by the continuity of V^{u), V~{u), V'^{u) and the 
compactness of {(m, r’*',r“, r°) G A x R'” x M^/maxdlr+H, ||r“||, ||r°||) = l}, there 
must be a (u, r~,r^) such that max(||r+||, ||r“||, ||r°||) = 1 and r+ + r“ + r° = 0 
which contradicts assumption (3.1). Thus: 



II 


< 


/3 



(3.3) 


The stable, unstable and neutral subspaces are also invariant under the differential 
of the map (ps{',h), i.e., if Ui+i = (ps{ui,h) and Vi+i = {Dips{ui,h))vi, then 


( V+im) v^+i G V+im+i), 

< 'Ci G V~{u^) Vi+i G V~{ui+i), (3.4) 

{ Vi € V°{ui) Vi+i G V°{ui+i) 

Because of their relative simplicity, studies of uniformly hyperbolic dynamical 

systems (also known as "ideal chaos") have provided a lot of insight into the proper¬ 
ties of chaotic dynamical systems [19]. Although most real-life dynamical systems are 
not uniformly hyperbolic, they can be classified as quasi-hyperbolic: results obtained 
on hyperbolic systems can often be generalized to them [20]. This proof covers the 
convergence os LSS for uniform hyperbolic flows, nevertheless, numerical results have 
shown that the algorithm also works for non-ideal chaos [14[. 


4. Neutral subspace and non-uniform discretization of trajectories. 

One should bear in mind that the dynamical system is continuous which means that 
a solution {u{t),t G to equation (1.1) forms a continuous trajectory in phase 
space. Thus, the sequence {ul^’°°\i G N} is no more than a sequence of sample 
points on the continuous trajectory, the time stepsize between two consecutive points 
being h. The neutral subspace V^{ui) which is unidimensional, is constituted by the 
line tangent to the continuous trajectory at the sampling point Ui. 

Consequently, a perturbation in the direction of the neutral subspace around the point 
Ui can be interpreted as a time shift. This means that for i G N* and e infinitesimal, 
if u{ti) = Ui : 


u{ti + €) = Ui + ef{ui,s) 


(4.1) 


and 


u{ti -I- e) = (ps{ui-i,h -I- e) = ips{ui-i + €f{ui-i,s),h) (4.2) 


which implies : 


(ps{ui-i + ef{ui-i,s), h) =Ui + ef{ui, s) 


i.e. 


{Dq:>s{ui-i,t))f{ui-i,s) = f{ui,s) 


(4.3) 


(4.4) 


Since ||/(ui,s)|[ < supj,gg(|[/(-, s)|[^) < oo for all i G N*, then (4.4) implies that a 
small perturbation in the direction of the neutral subspace remains bounded under 
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the action of forward or backward iterations. On the contrary, a small perturbation 
in the stable or unstable subspace gets amplified exponentially under the action of 
backward or forward iterations respectively. 


The second point to be discussed in this section is the fact that the discretization 
of a continuous trajectory doesn’t have to respect a uniform stepsizing. A better way 
to discretize a trajectory would be the following one : {{ui^Ti),i € N} is a sampling 
of a continuous trajectory if Ui = ^ps{ui-l^ hri) where G R. In order to have : 


{JM 


L h J 


lim 
T —)- + oo 

2=0 


^3 3 


0{h) 


(4.5) 


we only need sup(/iri) —>■ 0 as h —^ 0 which actually happens if, for instance, the time 

[Tl 

dilatation factors are bounded means ^^e above expression). From 

now on, a discretization of a trajectory is a sequence of couples (ui,Ti). 


5. Structural stability and the shadowing direction. In this section, we 
will prove a variant of the shadowing lemma for the purpose of defining the shadowing 
direction and prove its existence and uniqueness. 

The hyperbolic structure ensures the structural stability of the attractor A under 
perturbation in s [21], [22]. Without loss of generality, we will assume that s = 0 and 
choose a sequence G Z} such that : 

M°+i = >Po{u°,h) 

In this case, the superscript in is the value of the parameter s. h and T = oo are 
fixed so they do not appear in the notation. 

Theorem 5.1 (Shadowing trajectory). If the system is uniformly hyperbolic 
and ifs continuously differentiable with respect to s and u, then for all sequence 
{m°, z G Z} C a satisfying it? = (^o(w°_i, h), there is a M >0 such that for all jsj < M 
there is a sequence {(itf,Tf),i G Z} C R™ satisfying jjitf — it?|| < M, ||Tj®|| < M and 
itf = ^Ti) for all i G Z. Furthermore, itf and rf are i-uniformly continuously 

differentiable to s. 

The i-uniform continuous differentiability of itf means that for all s G {—M, M) and 
e > 0 there is a 5 > 0 such that if |s — s'j <5 then ||^)^(s) — -^(sOII < ^ 
I^W-$(^')l<eforalH. 

To prepare for the proof, let B be the space of bounded sequences in R”’' and Vi 
the hyperplane of R™ defined by Vi = F+(i(?)©F“(it?). We introduce V as the space 
of bounded sequences {vi,i G Z} such that Vi G Vi for all i G Z {vi has no components 
in the neutral subspace). In other words : 

v = ([[V)nB 

iei. 

Finally, by considering the space T of bounded sequences of R, we denote A the 
product of V by T : 


A = V X T 





7 


We then introduce the notation (v,x) = {{vi,Ti),i € Z} £ A where v £ V, r £ T 
and dehne the norm : 

||(v,r)||A = sup(||r;i||) + sup(|Ti|) = ||v||oo + ||r|loo 

zGZ zGZ 

As dehned above, the space A is a Banach space. 

We can now define the map F : A x R —>■ B as : 

V(v, t) £ a, Vs £ R, F{{v,t),s) = + Vi- hTi-i)),i £ Z} 

For a given s, F{{v, r), s) = 0 if and only if {(rt? +Vi), i £ Z} samples, with timesteps 
{hTi,i £ Z}, a continuous trajectory satisfying (1.1). We use the implicit function 
theorem to complete the proof, which requires F to be differentiable with respect to 
(v, x) and its derivative to be non-singular at v = 0, r = 1 and s = 0. 


Lemma 5.2. Under the conditions of theorem 5.1, F has Frechet derivative at all 
(v,x) £ A and |s| < M: 


(DF((v, t),s))(w, e) ={wi - {D(ps{u°_^ + Vi-i,hTi_i)){wi-i) 

- hei-idh^s{u^i-i + Vi-i,hTi-i),i £ Z} 


where (w, e) £ A 
Proof. We have: 

F((v wf, T e), s)) - j^((v, r), s) ^ f Wj 

ll(w,e)||A l||(w,e)||A 

V’s{Ui_i +Vi-i +w^-i,h{Ti-i +ei-i)) -(ps{u°_i +Vi-i,hTi_i) 


(w,e)||A->'0 


IKW,ej||A 


w,e 


hei-idh(psiu°_^ +v^-i,hn-i) 


w,e 


)} 


-I- Vi-l,hTi-i)){wi-l) 

ll(w,e)||A 


} 


(5.1) 

(5.2) 

(5.3) 

(5.4) 

(5.5) 

(5.6) 


in the A norm thanks to the uniform continuity of Dip and dhP on the compact set 
A. Now, we only need to prove that the linear map we obtained is bounded. Since 
Dp and dhP are continuous they are uniformly bounded in the compact set A. Thus, 
the norm of the linear map is less than (1 -|- ||Z1(/3||^ -|- H^/iV^H^). □ 


Lemma 5.3. Under conditions of theorem 5.1, the Frechet derivative of F at 
(v,x) = (0,1) and s = 0 is a bijection. 

Proof. The Frechet derivative of F at (v,t) = (0,1) and s = 0 in the direction 
(w, e) is : 

{DF{{0,l),Q)){w,e) = {wi - {Dpo{ui_j^,h))wi-i - he^-l^hPo{u^_J^,h),i £ Z} 


To prove its bijectivity, we only need to show that for any r = {ri} £ B there is a 
unique (w, e) £ A such that {DF{{0, 1), 0))(w, e) = r 

In this case, we can find an analytical expression for the pre-image of r. Let (w, e) 
be defined as : 




{^iDpoiu°_j, jh))r^_. - '^{Dpoiu°^j,-jh))r+^j, 
j=o j=i 


1 {r°^i;dhPo{u°,h)) \ 

h {dhpo{Ui, h); dhPoiu°, h))) 

(5.7) 


I 










We can verify that Wi — {D(Pq{u^_^, h)){wi-i) — hei-idh^o{u^_^, h) = for all 
We still have to ensure that (w, e) belongs to A. Based on (3.4), we notice that the 
Wi we’ve just defined belongs to 14 = V^iui) ® y~{Ui)- Then, thanks to (3.1), we 
can write = rf +r~ +r° where rf G l/+(u°), r~ G V~{u^) and G F°(m°). Since 
V'^{u)^ V~{u) and V^{u) are continuous with respect to u and A is compact: 

max(||r+||, ||r“||, ||r°||) < for all i (5.8) 


where /3 > 0. 
Consequently, for all i: 




i=o 

OO 


i=i 


i=o j=i 

2C IkllB 


< 


1 - A'* /3 


(5.9) 

(5.10) 

(5.11) 


because and V are invariant under Dipo (property (3.4)). Thus, Wi is uniformly 
bounded and w G V. In the same way, we show that for all i: 


. < l \\rU\\\9nMulh)\\ ^ |k°+i|| 

“ /i \\dhipoiUi,h)\\'^ ~ h\\dipo{_u°,h)\\ 

^ IkllB 

~ hPra 


(5.12) 

(5.13) 


where m = inf^gA {/("^i 0)} >0. Consequently, is uniformly bounded which leads 
to e G T and (w, e) G A. 


Because of linearity, uniqueness of (w, e) such that {DF{{0, 1), 0))(’w, e) = r only 
needs to be proved for r = 0 . Since R™ = lA+(u°) 0 V~{u^) 0 = 0 is 

equivalent to rf = r~ = = 0. Thanks to property (3.4), by splitting Wi = wf +w~ 

and knowing that hei-idh'^o{u^-i^ h) G y°(u°) , we have: 

0 = rt +r~ = {w+ - {D<fh{u°_i,0))w+_^) + (wk - {D<fh{ui_i,0j)w-_^) (5.14) 

where the two parentheses are in V'^{ui) and V~{u^) respectively. Again knowing 
that R™ = 0 V~{u^) 0 14°(u°), both parentheses should be equal to zero. 

This is true for all z, so we obtain : 

wf = {D^po{u°_^,h))w+_^ = ... = {Dip''Q~^\u°,h))w+ (5.15) 

u!~ = {D(fo{Ui_i, h))w-_^ = ... = h))w- (5.16) 

for all j < i. Based on the properties of uniform hyperbolicity, ||w^|| < 
and llzukll < CX^^^~^^\\wJ\\. If for some j we have w'j' k 0, then: 

> Ik+ll > ^ ^ Ik+ll for all i > j (5.17) 

^ Based on the fact that D(Pq{u^_^^ h)DipQ(u^_j_^, jh) = 

Dvo{Ui_i, h)Dipo{u9_^_^., -jh) = D(/Jo(«i_i+j' (-1 + 1)^) 


Dipo{u^_j_l,{j + l)h) and 
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which means that {wi^i G Z} is unbounded (0 < A < 1). Similarly, if Uj ^0 for some 
i then: 


lb /;-11 , 

—^ > ||w+|| > ^ ||w+|| for all j < i (5.18) 

and {wi,i € Z} is also unbounded. Consequently, for {wi} to be bounded we must 
have Wi = wf + w~ = 0 for all i. 

On the other hand, showing that = 0 is trivial: 

0 = r° = -hei-idhipoiu°, h) (5.19) 

Since Wdh^Poiu^, h)\\ > m > 0 then ei_i = 0. This is true for all i, which means that 
e = 0. This proves the uniqueness of (w, e) for r = 0. □ 


Proof, (of theorem 5.1) Since F((O,l),0) = = 0, (0,1) is 

a zero point of T at s = 0. Based on this information and on the two previous 
lemmas, the implicit function theorem states that there exist M > 0 such that for all 
|s| < M there is a unique (v'^jT®) satisfying ||(v®,t®)||a < M and T((v®,t®),s) = 0. 
Furthermore, this (v®,r®) is continuously differentiable to s, i.e., ^ G A is 


continuous with respect to s in the A norm. By the definition of derivatives (in A), 


d(v‘ 


ds 


= Continuity of 


d(v®,T®) 
ds 


in A then implies that and 


are i-uniformly countinuous with respect to s. By defining: 


{«,T®),i G Z} = {(w° + : 


’),*GZ} 


(5.20) 


we finally obtain the results of theorem (5.1). □ 


If we return to the expanded notation of uf, this theorem states that for a dis¬ 
cretization 1)} of a continuous trajectory satisfying (1.1) for s = 0, there 

is a series also being a discretization of a continuous trajectory 

at nearby values of s. In addition, shadows I)} mean¬ 
ing that r/’^^’°°^)} is close to 1)} when s is close to 0. Also, 

^^s{h,oo} ^^s{h,oo} 

{( —)} exists and is ^uniformly bounded. 


The shadowing direction is defined as the uniformly bounded 


series : 


{( 


ih,oc)\ { h,oo\\ 1 

v} ' ,77^ 


))-{( 


du. 


s{/i,oo} 


ds 


dr. 


s{/i,oo} 


s=0 


ds 


s=0 


)} 


G A 


(5.21) 


In addition, we can find 2 constants and independant of h such that 

for all i\ 


^ ||v{oo}|| < ||J7^°“>|| (5.22) 


Please refer to appendix A to see how these constants are built. 
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6. A simpler result. In this section, we prove an easier version of Theorem 
LSS in which we replace the solution to the con¬ 

strained least squares problem (2.3) by the shadowing direction we found earlier 


^ ...JT]} which can be written 

h is to be displayed explicitly. 


I 'A 


):* = if 


Theorem 6.1. If uniform hyperholicity holds and Lps{-,h) is continuously differ¬ 
entiable for all h, then for all constinuously differentiable function J : R™ X R —^ R 
whose infinite time average : 


{J){s) 


lim 
T ^ + oo 


1 

T 



J{u{t), s)dt 


where 


du 

dt 


f{u,s) and u{0) = uq ( 6 . 1 ) 


is independant of the initial state uq, let {(yj^’°°\r]j^’°°^),i = 1 ,...,[^]} be the se¬ 
quence of shadowing direction, then: 


[T] 

+ {dsJ{ui,Q)) + (?7f^’'^f(J(M„0) - (J)(0))) 

i—1 

( 6 . 2 ) 

, [?1 

^ T^h^or'^ + idsJ{Ui,0)) -I {iqj^’°°'^{J{u„0) - (J)(0))) 

(6.3) I 


Proof. The proof is essentially an exchange of limits through uniform convergence. 
Since (J) is independant of uq, we set uq = as defined in the previous section 

and we know that hT^^^'°°^). In order to simplify the notations, 

^s{/t.oo} ^s{h,co} expressed as rtf and rf in what follows. The integral could 

be discretized into an infinite sommation that involves the time dilatation factors: 




lim 

T—f-\-oo 


1 

T 


J{u{t),s) 


lim lim > 

h^O T->+oo 


r!J{ut,s) 

3 



lim lim y 
T—f-\-oo h—¥0 ‘ ^ 
i—1 


TfJ{uf,s) 
^3 3 


(6.4) 

(6.5) 
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The two limits can be permuted at will but we will only keep one notation for the 
remaining of the proof. We can write: 


d{J) 


ds 


s=0 


= lim 


(J)(.)-(J)(0) 


= limlim lim T W'■n' 4 «l,») 

.<?—A>—i.n T—i.-1-rv-i \ G f ^ 


s-s-0/t-j-O T-H-oo \ S —' h^-T, 

' 7=1 


1 3 


T 


= lim lim lim , / , 

s —^0 h —^0 T —^+oo \ 5 h/ 

2 = 1 




'3 '3 ^3 3 ^3 '3 




T 


= lim lim lim 

s^OT—>-+oo \ 5 


1 ^ h{J{uls)- J{ul 0) + (r/ - 1) J«, s)) 


2 = 1 




3 3 


- 1)) X W(m°,0) 


T{hY.,r^ 


= lim lim lim 
s —^0 h —^0 T —^-j-cxD V 5 


1 L h J 

7E 


h{J{uf,s) - J(m°, 0) + (r/ - l)J(Mf,s)) 


T + /.E,(r/-l) 


- 1)) X W(m°,0) 
nr+hj:Ar^ -1)) 


= lim lim lim 

s —^0 h —^0 T —^-j-OO y 2”* 


^ s) — J{u°,0)) ^ 


s-^0 h->-0T^+oo \ IS 1 ) 


Let us eliminate lims_>o in the first term. We define : 

ll = ( 6 . 6 ) 

Then, thanks to the mean value theorem, for all i there exist an ^i{s) G [0,s] such 
that: 


(J«,s) - J{u°,0)) ^ 
s 

Consequently: 


(6.7) 


lim lim lim 
s—)-0 T—>-+oo 



iJ{uf,s) - J(m°,0)) ^ 


lim lim lim . / 

s—>^0 h—^OT ^+oo \ T 

2=1 




(«) 


( 6 . 8 ) 

We can choose a neighborhood of Ax {0} that contains (uf, s) for all i (for s sufficiently 
small) and in which both {DJ(u,s)) and 9sJ(u, s) are uniformly continuous. Since 
the are i-uniformly continuous and bounded, for all e > 0 there exists M > 0 
such that for all |^| < M: 


Il 7 z— 7 *°ll<e Vi 




















12 


Thus, for all |s| < M, |^i(s)| < \s\ < M for all i, therefore for all h, T {h << T) : 


h 

T 



i=l 







i=l 


7°ll < 


(6.9) 


Hence, 


lim lim 

h —^0 T —^-1-C30 



— lim lim 
h —^0 T —>^-l-C30 



< e 


( 6 . 10 ) 


Finally, 


lim lim lim 

s —^0 h —^0 T —^-|-oo 



lim lim 

h—¥0 T —>-+co 



( 6 . 11 ) 


which grants us the desired result for the first term via the definition of 7 °. 


For the second term, J is continuously differentiable thus continuous and the 
(wf,r®) are i-uniformly continuously differentiable and bounded. Based on that, for 
s sufficiently small, we can find a compact neighborhood of A x {0} that contains 
(uf,s) for all i G Z and in which J{u,s) will be uniformly continuous. Consequently, 

the sequence { J{uf, s), i G Z} which can be writen { K j{uf,s),i S Z} 

^ {fc.CXj} 

converges uniformly to — J{u^,s),i G Zj when s goes to 0. Because the term 
T 0 ) does not depend on s at all, we finally have: 


lim lim lim 

s—>-0 /i—>-0 T—>-+oo 


Krt - 1 ) 

Ts 


{j{ul,s)--^hJ{u°,Q)) 


lim lim 

h —^0 T —^-j-oo 


fhrl^ 


(J(«°,0)-(J)(0)) 


which conludes the proof. 

□ 


7. Computational approximation of the shadowing direction. The main 
task of this section is to provide a bound for : 


{/i,oo} 

el = Vi — vl 

{h,T} {h,T} 

e) = W - W 


(7.1) 

(7.2) 


for i = 1 ,..., ^ where the 


are the solution to the least squares problem: 


i=l 

s.t. = {D(ps{ui, h))vj^’'^^ + {ds(ps{ui, h)) + dh(psiui, h), 


(7.3) 

(7.4) 


The shadowing lemma guarantees the existence of a shadowing trajectory, but pro¬ 
vides no clear way to compute This section suggests that the 

solution to the least squares problem gives a useful approximation of the shadowing 
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trajectory allowing us to compute^^^. Without loss of generality, we consider that 
s = 0 in (7.4). By definition, the shadowing trajectory satisfies: 

) (7.5) 

After taking the derivative to s on both sides for s = 0, we obtain: 

^ h))4^’°°^ + + h7^|^’°°^^h‘fio(u^, h) (7.6) 

Thus, the shadowing direction satisfies the constraint in equation (7.3) and: 


(7.7) I 


Consequently : 


max < 


• {h,T} {h,T} {/i.oo} 

Since = Vi ~ Vi ■ 




(7.8) 


max II11 < 


r/||v{°°>||| 


(^ 


+ ||^{-}||2j%||^{oo}| 


(7.9) 


We can get a similar bound for max, ||ej^^’^^°|| : 


m{?..T}0m/ /r(||vi°°}||2 +a||,y{-}||2)i ||v{->||B 

max e, < ‘'-^- 


(7.10) 


Concerning max^ ||e|^’^^~''|| and max^ ||e^^’^^ ||, combining the constraint equa¬ 
tion (7.4) as well as (7.6) we obtain : 


\ epp~ = Dipo{ui,h)el^’'^^~ 

Thus, since G V'^{ui), G V~{ui) and based on the following 

definition of the unstable and stable subspaces : 


y+(u) := {u G R™/ ||(£»vj,(M,t))u|| < C'A-‘||?;||,Vt < 0} (7.11) 

t/-(u) := {z; e R™/ ||(£)(^,(M,t))w|| < C'A‘||u||,Vt > 0} (7.12) 

we obtain : 

max lle^^’^^^ll < max(l, C')||ef^;^^^|| (7-13) 

i IttJ 

max ||e|^’^^~|| < max(l, (7)116^’^^” || (7-14) 
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In addition to that, (/5 o(m, 0) = u which means that Dlpq{u^ 0) = / where / is the 
identity operator in R™. For u G A, h G H and s G S where H and S are compact 
sets containing h = 0 and s = 0, we can find a positive constant K such that: 

\\dhDipoiu,h)\\<K (7.15) 

since (•,•) is C^. dhD(pQ{u,h) belongs to a space of finite dimension so all norms 
are equivalent and any choice of norm || • || is valid. For a sufficienty small h such that 
hK < 1, we get : 


> (1 - hK)\\v\\ 

for all V G R™. Consequently : 


I h j I fc i 

Jh,T}-n 2 




2(i-l)nJh,T}-n2 


> l|ei 


{h,T}-.,2 ^ ~ (1 ~ hK)^^T:'^ 
" h{2K-hK^) 


which means that: 




< 


El 


Eii4' 


h{2K — hK‘^) ^ ^ ^ IIII 2 


1 — e 


-2TK 


1 — e 


-2TK 


EII4^ 


For T sufficiently big, e ^ which means that: 


||er'-|p<4M-y:i|el‘U-||2 


2=1 

Since then: 

Ilefu-I,. < 2(||„(‘U-||2 + < 2(||„i>..n||. + ||„f.=<.>||2) 

leading to: 




<^xE2|M->|||+a|M->|p) 

7 h. 


Finally : 


max lie,- ' " || < x max(l, C) x ^- x (2 ||v^°“^||b + a||r 7 ^°°^|p)y 


Ah,T}-. 


< eVt 


(7.16) 

(7.17) 

(7.18) 

(7.19) 

(7.20) 

(7.21) 

(7.22) 

(7.23) 

(7.24) 

(7.25) 

(7.26) 
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where E = max(l, C) x x (2||v'f°“^||| + ^ is a constant that doesn’t 

depend on h nor on T. In the same way we obtain : 

max||e|'‘’^^+|| < (7.27) 


Even though the bounds we found for and may 

depend on T and/or h, we will see in the next section that they are strong enough to 
prove the convergence of the algorithm. Furthermore, experimental simulations have 
shown that and doesn’t necessarily grow when ft, —>■ 0, T —oo and stay 

bounded [18]. 

8. Convergence of least squares shadowing. In this section, we use the 
results obtained previously to prove our initial theorem : 

Theorem 8.1 (THEOREM LSS). For a map V? (-,-) and a cost function 
J, the following limit exists and is equal to: 


d{J) 

ds 


lim lim 

/i—>^0 T—^oo 




+ (5sJ(Mi,s)) + (r 7 ,^^’^^(J(ui,s) - (j)(s))) 


lim lim — 
T ^oo T 




+ {dsJiui,s)) + s) - (J)(s))) 


I 


Proof. Because J is and A is compact, {DJ{ui, 0)) is uniformly bounded, i.e., 
there exists a constant A such that ||Z?J(ui, 0))|| < A for all i. Let be defined 

as in the previous section, then: 


tYI [iddJ{ui,s))vj'"’'^^ + {dsJ{ui,s)) + - (J)(s))) 

h r 

~tY [(ddJ{ui,s))vl^'°°^ + {dsJ{ui,s)) + {r]l^'°°'^{J{ui,s) - (J)(s))) 

2=1 


( 8 . 1 ) 

( 8 . 2 ) 

(8.3) 


7 LTTJ 

fl 

2=1 




{h,T}+ Jh,T}- 




(8.4) 


< 




2=1 


fY [iddJ{ui,s))ej^''^'^° + (e,^^’^^(J(wi,s; 

(8.5) 


2=1 


- {J)B))) 


For the first term: 
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T 




i=l 


<- 


T 


^||pJK,s))||||ef’^>+|| 


Z = 1 


-^\\iDJ{u.,. 


Ah,T}-: 


( 8 . 6 ) 

(8.7) 




< 


i=l 

h 2AC 
r(i-A'*) 




i;Vr 


1 2y4C£; 

-7f 


( 8 . 8 ) 

(8.9) 

( 8 . 10 ) 


which goes to 0 when T increases. Thus, we notice that, in the stable and unstable 
subspaces, the difference between the shadowing trajectory and its ap¬ 
proximation decreases extremely fast so that the whole term | y X)l=i {DJ{ui, s))(ej-^’^^~''-|- 

e|^’^^~)| tends to 0. 


The situation is more complicated for the second term since there is no reason 
for and to decrease when ^ increases. The cancellation of the second 

term is the result of the mutual cancellation of the elements in the sommation as 
we shall see. Based on the sampling points t)*^^’°°^)} for the continuous 

shadowing trajectory found in section 5, we consider the new set of sampling points 

which satisfy the following relation : 


lim 

s—»-0 



_ AAT} 


( 8 . 11 ) 


for all i. We can notice that the new sampling points describe the same continuous 
trajectory as the old set of values. Assuming that and are bounded, we 

have for a sufficiently small s: 


ITTI 

(J)(s) = lim lim 

T-s-l-oo ^ 
i=0 


hjr! + sef'^^)J{uf,s) 


( 8 . 12 ) 


We would obtain by following the same operations we did in section 6 (but upside 
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down this time) : 


lim lim s) - (J)(s))) 

h^OT^oo \ 1 ^—' L 


= lim lim lim ( — > 
h^O T-S -+00 s-)-0 T 

'he} 




i=0 

{h,T} 


+ a rST» I*?. (-\- - f E 


= lim lim + =) 

)-0 T^+oo s—>-0 V 5 


2=0 






) 


1 + s4^’^^)>^(w?,s) hTfJ{uf,s) 


= lim lim lim ( — 

s —V n h _^ n '"P _^ -L ("v"! \ c ' ^ ^ ^ 

+ 0(s) =0 


.^o^^oT^+oc Vs^ +sef’^}) 

(J)(s)-(J)(s) 


= lim 

S —)-0 


0{s) 

0{s) 


(8.13) 

(8.14) 

(8.15) 

(8.16) 

(8.17) 

(8.18) 


However, this is not necessarily true since and may grow as ^ 

increases. Permuting the limits is much more delicate but is still possible. Please 
refer to appendix B for further details about how to permute the limits. The idea is 
to use the relation : 




2=1 


2=1 


|2 

Ib 


(8.19) 


) 

I 


to show that most remain bounded and that the contribution of the 

unbounded terms fades out. 


In conclusion : 


d{J) 

ds 


y E + {dsJ{Ui,s)) + - (J)(s))) 

2=1 

( 8 . 20 ) 


d{J) 

ds 


, [xl r 

y E [(^d{ui,s))vl^'°°^ + {dsJ{ui,s)) + - (J)(s))) 

( 8 . 21 ) 


+ 


h r 

^E + {^sJ{u^,s)) + (? 7 ,^^’^^(J(u„s) - (J)(s))) 


h 

f 


m ^ 

E [{DJ{ui,s))v^^'°°^ 


+ {dsJiUi,s)) + (77,^^’°°^(J('Ui,s) - (J)(s))) 


( 8 . 22 ) 

(8.23) 

(8.24) 
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and we have showed that both terms go to 0 as T —>■ oo and h —>■ 0. This concludes 
the proof. □ 


9. Conclusion. As we have shown through this paper, LSS gives us a good 
estimation for when the dynamical system is uniformly hyperbolic. After running 
a simulation for a given s, an arbitrary initial condition uq and an uniform time 
steptize of h, we obtain a sequence of reference sampling points^ {u|, z = 1, • • • , ^ 

If we had access to the shadowing direction, we would easily compute : 


d{J) 

ds 


h 

T 


m 


(9.1) 


However, in real-life problems we usually do not have access to the stable and unstable 
subspaces around each uf prohibiting the usage of the closed form expression of 
and Thus, we have no other choice than computing an approximation of the 

shadowing direction. This approximation is given by the solution to the least squares 
problem: 


s-t. vj'l’P = {D(ps{ut, h))vj^’'^^ + dsg}s{ul, h) + dhifsiK, h), 


(9.2) 


After solving this quadratic optimization problem, we estimate using expression 
(9.1) again where the 77 !^’°°^) are replaced by 77 ^^’^^). As we have 

seen previously, this estimation converges to the real value of when the time 
stepsize h is refined and the integration lapse T increases. 
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Appendix A. In our case, we have an explicit expression for {v} 


i 


oo 


oo 



(A.l) 


{h.oo} ^ 1 {ds(poiu°,h)°;dh^oiu°,h)) 

h {dh(poiu°,h)]dh(poiu°,h)) 


(A.2) 


As we did previously: 


CXD OO 



(A.3) 


OO 


OO 



(A.4) 


3=0 

2C sup„gA(||ds(po(M,/r)||) 

1 - A'' p 


< 


(A.5) 


I 
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so : 


{/i.oo},, ^ 2C sup„gA(ll^«^o(»)^)ll) ^ 2C sup^^j^{\\ds>po{u,h)\\) 

P “ /31n(f) h 

^ 2C f^^ds(po{u,h) - ds(fo{u,Q) 

0(^ 

- R] (<=\ {\\dhdsMu,h)\\) = ||v^°°>|| 

(u,h)^KxH 


(A.6) 

(A.7) 

(A.8) 


where 77 is a compact set of R’*' (for example [0,1]) and ds<fo{u^ 0) = 0. We’ve also 
used the taylor expansion of 1 — for /i —>■ 0 and assumed that duds is well-defined 
and continuous on the compact set A x H. 


Following similar steps, we obtain : 

{/I.oo} ^ SUp„gA(||5s<Po(M,/t)ll) 

^ “ Phm 

^ suP(«,;»)gAxg {\\dhds<fo{u, fe)||) 

“ /3m 


\V 


{oo}| 


(A.9) 

(A.IO) 


Appendix B. Let C G K"*" be an arbitrary bound and we will assume that a = 1 
for simplicity reasons (without loss of generality). If Ue is the number of elements 
that are bigger or equal to C, we have: 


He < 


/lC2 


(B.l) 


with the equality being verified in the worst case scenario where all the unbounded 
terms are equal to C, all the bounded terms are equal to 0 and UeC^ is exactly equal 
to ^(||vt°°5^ llg -I- 11??’^°°^IP). Then, let us compute the contribution of the terms that 
are unbounded. For that purpose, we introduce an indicator function S that is equal 
to 0 when both and are below C and equal to 1 when at least one of 

them is bigger (or equal) than C. We have: 


lim lim (-^A(i)[(7AJ(ui,s))e|^’^^° + (e|'*’^^(J(Mi,s) - (J)(s))) 
h^O T—)-oo \ 1 ^' 


2=1 




h—^OT—^oo \ T 


7||A 


.)) 


2=1 


[-] 

- (^max(||7AJ||^,2||J||^) ^(5(i)(||ef’^^°|| + 

- f|;max(||7AJ||^,2|lJ||^) x UeC 

?i—s-O T—t-oo y J 

- f|:max(||7AJ|l^,2|lJ||^) x ^ 

h—^O T-700 A 7 

^ max(||7AJ||^,2||J||^)(||v{-}||| + ||^{-}|p) 

C 



(B.2) 

(B.3) 

(B.4) 

(B.5) 

(B.6) 

(B.7) 
















21 


We have used the fact that X]l=i i® maximized when we 

have the maximum number of unbounded elements and when all of them have the 
same value. The worst case scenario is again the one where we have Ue unbounded 
elements all equal to C. Since C is arbitrarly big, the contribution of the unbounded 
terms is as small as we want. Then, the limits can permuted in the expression : 


lim lim 
/i—>-0 T—^oo 


2=1 





{J{uf,s) - J(Mf,s)) ^ 


+ lim lim lim 

h —^0 T —^-j-cjo s —^0 


= lim lim lim 

h —^0 T —^-|-oo s —^0 


= lim lim lim 

s—>0 h—^0 T—y+oo 


Ah,T} 




2=0 


( 1 -«(<))-!: 


2 = 0 


3 

.S))j 

+ se|^’^^)J(u'*,s) 

hTiJ{uf,s) 



Kt! + se\^''^^)J{Ui^,s) 

hT^J{uf,s) 




(B.8) 

(B.9) 

(B.IO) 

+ 0(s) 

(B.ll) 


s->0 S 


+ 0(s) = 0 


The last equality comes from property (4.5) we had on Riemann sums. 


+ 0(s) 
(B.12) 
(B.13) 











